##############################################################################################################################
## Fatality Thresholds, Causal Heterogeneity, and Civil War Research: Reconsidering the Link Between Narcotics and Conflict ##
##############################################################################################################################
################################################ Noel Anderson & Alec Worsnop ################################################
##############################################################################################################################
###################################################### Code for Tables #######################################################
##############################################################################################################################

## Install required packges
install.packages(c("foreign","cem","mvtnorm","MASS","lmtest","MatchIt","optmatch",
                 "sandwich","lme4","car","fields","stargazer","survival"))

## Load libraries
library(foreign)
library(cem)
library(mvtnorm)
library(MASS)
library(lmtest)
library(MatchIt)
library(optmatch)
library(sandwich)
library(lme4)
library(car)
library(fields)
library(stargazer)
library(survival)

## Load the dataset
dataset <- read.dta("<<insert file location>>/dataset.dta")

## Create function to cluster standard errors
cl <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) 
}

###########################################
##### APPENDIX C (Regression Results) #####
###########################################

## Model 1: Including all observations
# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","ccode","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- dataset[keep.vars]
datafull <- na.omit(all.data)

## Run the model
modelfull=lm(BTTLDtotLN ~ DRUGP + ALLGEMP + hydrocarbonP + OUTZhydrocarbonP + 
               mountain + confbord + Weakrebel + Equalstrenght + coldwar + 
               internatz + popLN + ethfrac + ethPOL + DEM,data=datafull)

# Pull out clustered robust t-scores 
tvaluefull=abs(round(cl(datafull, modelfull, datafull$ccode), 2))[,3]

# Pull out clustered robust p-values 
pvaluefull=cl(datafull, modelfull, datafull$ccode)[,4]

# Exponentiate the coefficients
modelfull$coefficients=round(exp(modelfull$coefficients),3)

## Model 2: Excluding Conflicts with Drugs and <50 Battle Deaths
# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","ccode","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Remove 4 observations with drugs and less than 50 battle deaths
datatrim=data[!(exp(data$BTTLDtotLN)<50&data$DRUGP==1),]

# Run the model
model=lm(BTTLDtotLN ~ DRUGP + ALLGEMP + hydrocarbonP + OUTZhydrocarbonP + 
           mountain + confbord + Weakrebel + Equalstrenght + coldwar + 
           internatz + popLN + ethfrac + ethPOL + DEM,data=datatrim)

# Pull out clustered robust t-scores 
tvalue=abs(round(cl(datatrim, model, datatrim$ccode), 2))[,3]

# Pull out clustered robust p-values 
pvalue=cl(datatrim, model, datatrim$ccode)[,4]

# Exponentiate the coefficients
model$coefficients=round(exp(model$coefficients),3)

# Make table
rownames=c("Drug Cultivation (In zone)","Gem Production(In zone)","Hydrocarbon Prod. (In Zone)","Hydrocarbon Prod. (Out Zone)","Mountainous terrain (log)","Conflict at border","Weak rebel group","Equal strength","Cold War","Internationalized conflict","Population size (Log)","Ethnic fractionalization","Ethnic polarization","Democracy")
stargazer(modelfull,model, 
          se=list(tvaluefull,tvalue),
          p=list(pvaluefull,pvalue), 
          type="text",report = "vc*s",digits=2,
          model.numbers=FALSE,
          omit = "Constant",omit.stat = c("adj.rsq","ser","f"),
          covariate.labels=rownames,
          dep.var.caption="",
          dep.var.labels.include = FALSE,no.space=TRUE,
          column.labels = c("All Observations", "Excluding Four"))

###########################################
##### APPENDIX D (Regression Results) #####
###########################################

## Table D.1 
# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","ccode","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Generate thresholds and objects for each threshold in which observations are dropped 
btldt=seq(25,1000,25)
n.est=c()
t.est=c()

# Get the thresholds that drop observations
for (i in btldt){
  datatest=data[exp(data$BTTLDtotLN)>i,]
  n.est=c(n.est,dim(datatest)[1])
  t.est=c(t.est,sum(datatest$DRUGP))
}
drops=cbind(btldt,n.est,t.est)
drops=rbind(1,drops)
x=matrix(,dim(drops)[1],3)
for (i in 2:dim(drops)[1]){
  x[i,1:3]=ifelse(drops[i,2]==drops[i-1,2],NA,i)
}
drops2=na.omit(x[,1])
uniquen=(drops[drops2,])
uniquen=cbind(uniquen,1:29)

# Now run the regressions only for the new thresholds
for (i in uniquen[,1]){
  datatest=data[exp(data$BTTLDtotLN)>i,]
  formula = as.formula(paste(paste("BTTLDtotLN", '~'), paste("DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM",collapse='  ' ,sep="+")))
  model=lm(formula,data=datatest)
  # Pull out clustered robust t-scores 
  tvalue=abs(round(cl(datatest, model, datatest$ccode), 2))[,3]
  # Pull out clustered robust p-values 
  pvalue=cl(datatest, model, datatest$ccode)[,4]
  # Exponentiate the coefficients
  model$coefficients=round(exp(model$coefficients),3)
  # Number models, t values and p values.
  assign(paste("M",uniquen[uniquen[,1]==i,4],sep=""), model)
  assign(paste("T",uniquen[uniquen[,1]==i,4],sep=""), tvalue)
  assign(paste("P",uniquen[uniquen[,1]==i,4],sep=""), pvalue)
}

# Print results displayed in Appendix D.1
rownames=c("Drug Cultivation (In zone)","Gem Production(In zone)","Hydrocarbon Prod. (In Zone)","Hydrocarbon Prod. (Out Zone)",
           "Mountainous terrain (log)","Conflict at border","Weak rebel group","Equal strength","Cold War","Internationalized conflict",
           "Population size (Log)","Ethnic fractionalization","Ethnic polarization","Democracy")
colnames=as.character(uniquen[,1])

# First page
stargazer(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13, M14, 
          se=list(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14),
          p=list(P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14), 
          type="text",report = "vc*s",digits=2,column.labels=colnames[1:14],
          model.numbers=FALSE,omit = "Constant",
          omit.stat = c("adj.rsq","ser","f"),covariate.labels=rownames,
          dep.var.caption="Battle Death Threshold",dep.var.labels.include = FALSE,no.space=TRUE)

# Second page
stargazer(M15,M16,M17,M18,M19,M20,M21,M22,M23,M24,M25,M26,M27,M28,M29,
          se=list(T15,T16,T17,T18,T19,T20,T21,T22,T23,T24,T25,T26,T27,T28,T29),
          p=list(P15,P16,P17,P18,P19,P20,P21,P22,P23,P24,P25,P26,P27,P28,P29), 
          type="text",report = "vc*s",digits=2,column.labels=colnames[15:29],
          model.numbers=FALSE,omit = "Constant",
          omit.stat = c("adj.rsq","ser","f"),covariate.labels=rownames,
          dep.var.caption="Battle Death Threshold",dep.var.labels.include = FALSE,no.space=TRUE)

## Table D.2
# Create resources dummy for conflicts with hydrocarbons in and outside the conflict zone and with gem production
resources=c()
for (i in 1:dim(dataset)[1]){
  resourcesvar=ifelse(dataset[i,"OUTZhydrocarbonP"]==1|dataset[i,"hydrocarbonP"]==1|dataset[i,"ALLGEMP"]==1,1,0)
  resources=c(resources,resourcesvar)
}
dataset$resources=resources

# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","BTTLDtot","resources","ongoingCONF","Duration","DRUGP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","DEM","ccode","terr","gdplLN")  
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Create censored observations
data$ongoingCONF=1-data$ongoingCONF

# Run analysis
btldt=seq(25,1000,25)
for (i in uniquen[,1]){
  # Create dataset for each new battle death threshold
  datatest=data[exp(data$BTTLDtotLN)>i,]
  # Run model
  model=coxph(Surv(Duration,ongoingCONF)~DRUGP+resources+mountain+gdplLN+Weakrebel+Equalstrenght+internatz+
                coldwar+popLN+ethfrac+DEM+terr+confbord+cluster(ccode),data=datatest)
  # Number models
  assign(paste("M",uniquen[uniquen[,1]==i,4],sep=""), model)
}

# Print results displayed in Appendix D.2
rownames=c("Drug Cultivation (In zone)","Other Natural Resources","Mountainous terrain (log)","GDP/capita (Log)","Weak rebel group",
           "Equal strength","Internationalized conflict","Cold War","Population size (Log)","Ethnic fractionalization","Democracy",
           "Territorial Conflict","Conflict at border")

# First page
stargazer(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13,M14, 
          report = "vc*s",digits=2,column.labels=colnames[1:14],model.numbers=FALSE,
          omit = "Constant",omit.stat = c("adj.rsq","ser","f"),covariate.labels=rownames,
          dep.var.caption="Battle Death Threshold",dep.var.labels.include = FALSE,
          no.space=TRUE,type="text")

# Second page
stargazer(M15,M16,M17,M18,M19,M20,M21,M22,M23,M24,M25,M26,M27,M28,M29, 
          report = "vc*s",digits=2,column.labels=colnames[15:29],model.numbers=FALSE,
          omit = "Constant",omit.stat = c("adj.rsq","ser","f"),covariate.labels=rownames,
          dep.var.caption="Battle Death Threshold",dep.var.labels.include = FALSE,
          no.space=TRUE,type="text")

## Table D.3 - Intensity with no duration control (Table VI)
# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","popLN","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","BTTLDtotLN","ethPOL","DEM", "ccode")  
data <- na.omit(dataset[keep.vars])

# Run analysis
for (i in uniquen[,1]){
  datatest=data[exp(data$BTTLDtotLN)>i,]
  formula = as.formula(paste(paste("BTTLDintLN", '~'), 
                             paste("DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain",
                                   "confbord","Weakrebel","Equalstrenght","coldwar","internatz",
                                   "popLN","ethfrac","ethPOL","DEM",collapse='  ' ,sep="+")))
  model=lm(formula,data=datatest)
  # Pull out clustered robust t-scores 
  tvalue=abs(round(cl(datatest, model, datatest$ccode), 2))[,3]
  # Pull out clustered robust p-values 
  pvalue=cl(datatest, model, datatest$ccode)[,4]
  # Exponentiate the coefficients
  model$coefficients=round(exp(model$coefficients),3)
  # Number models, t values, and p values
  assign(paste("M",uniquen[uniquen[,1]==i,4],sep=""), model)
  assign(paste("T",uniquen[uniquen[,1]==i,4],sep=""), tvalue)
  assign(paste("P",uniquen[uniquen[,1]==i,4],sep=""), pvalue)
}

# Print results displayed in Appendix D.3 (Table VI)
rownames=c("Drug Cultivation (In zone)","Gem Production(In zone)","Hydrocarbon Prod. (In Zone)",
           "Hydrocarbon Prod. (Out Zone)","Mountainous terrain (log)","Conflict at border",
           "Weak rebel group","Equal strength","Cold War","Internationalized conflict",
           "Population size (Log)","Ethnic fractionalization","Ethnic polarization","Democracy")
colnames=as.character(uniquen[,1])

# First page
stargazer(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13, M14, 
          se=list(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14),
          p=list(P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14), 
          type="text",report = "vc*s",digits=2,column.labels=colnames[1:14],
          model.numbers=FALSE,omit = "Constant",
          omit.stat = c("adj.rsq","ser","f"),covariate.labels=rownames,
          dep.var.caption="Battle Death Threshold",dep.var.labels.include = FALSE,no.space=TRUE)

# Second page
stargazer(M15,M16,M17,M18,M19,M20,M21,M22,M23,M24,M25,M26,M27,M28,M29,
          se=list(T15,T16,T17,T18,T19,T20,T21,T22,T23,T24,T25,T26,T27,T28,T29),
          p=list(P15,P16,P17,P18,P19,P20,P21,P22,P23,P24,P25,P26,P27,P28,P29), 
          type="text",report = "vc*s",digits=2,column.labels=colnames[15:29],
          model.numbers=FALSE,omit = "Constant",
          omit.stat = c("adj.rsq","ser","f"),covariate.labels=rownames,
          dep.var.caption="Battle Death Threshold",dep.var.labels.include = FALSE,no.space=TRUE)

## Table D.3 - Intensity with duration control (Table VII)
# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","popLN","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","BTTLDtotLN","ethPOL","DEM", "ccode","durationLN")  
data <- na.omit(dataset[keep.vars])

# Run analysis
for (i in uniquen[,1]){
  datatest=data[exp(data$BTTLDtotLN)>i,]
  formula = as.formula(paste(paste("BTTLDintLN", '~'), 
                             paste("DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain",
                                   "confbord","Weakrebel","Equalstrenght","coldwar","internatz",
                                   "popLN","ethfrac","ethPOL","DEM","durationLN",collapse='  ' ,sep="+")))
  model=lm(formula,data=datatest)
  # Pull out clustered robust t-scores 
  tvalue=abs(round(cl(datatest, model, datatest$ccode), 2))[,3]
  # Pull out clustered robust p-values 
  pvalue=cl(datatest, model, datatest$ccode)[,4]
  # Exponentiate the coefficients
  model$coefficients=round(exp(model$coefficients),3)
  # Number models, t values, and p values
  assign(paste("M",uniquen[uniquen[,1]==i,4],sep=""), model)
  assign(paste("T",uniquen[uniquen[,1]==i,4],sep=""), tvalue)
  assign(paste("P",uniquen[uniquen[,1]==i,4],sep=""), pvalue)
}

# Print results displayed in Appendix D.3 (Table VII)
rownames=c("Drug Cultivation (In zone)","Gem Production(In zone)","Hydrocarbon Prod. (In Zone)",
           "Hydrocarbon Prod. (Out Zone)","Mountainous terrain (log)","Conflict at border",
           "Weak rebel group","Equal strength","Cold War","Internationalized conflict",
           "Population size (Log)","Ethnic fractionalization","Ethnic polarization","Democracy",
           "Duration (log)")

# First page
stargazer(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13, M14, 
          se=list(T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14),
          p=list(P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14), 
          type="text",report = "vc*s",digits=2,column.labels=colnames[1:14],
          model.numbers=FALSE,omit = "Constant",
          omit.stat = c("adj.rsq","ser","f"),covariate.labels=rownames,
          dep.var.caption="Battle Death Threshold",dep.var.labels.include = FALSE,no.space=TRUE)

# Second page
stargazer(M15,M16,M17,M18,M19,M20,M21,M22,M23,M24,M25,M26,M27,M28,M29,
          se=list(T15,T16,T17,T18,T19,T20,T21,T22,T23,T24,T25,T26,T27,T28,T29),
          p=list(P15,P16,P17,P18,P19,P20,P21,P22,P23,P24,P25,P26,P27,P28,P29), 
          type="text",report = "vc*s",digits=2,column.labels=colnames[15:29],
          model.numbers=FALSE,omit = "Constant",
          omit.stat = c("adj.rsq","ser","f"),covariate.labels=rownames,
          dep.var.caption="Battle Death Threshold",dep.var.labels.include = FALSE,no.space=TRUE)

#################################
##### APPENDIX F (Matching) #####
#################################

## Total battle deaths balance statistics (top panel, Table VIII)
# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","BTTLDtot","ccode","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- dataset[keep.vars]
data <- na.omit(all.data)

# Run analysis
matchrun=c()
model <- matchit(DRUGP~ALLGEMP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+
                     Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+
                     DEM, data=data, method="nearest") 
matchrun=cbind(summary(model)$sum.all[,c(1:2,4)],summary(model)$sum.matched[,c(1:2,4)])

model <- matchit(DRUGP~ALLGEMP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+
                     Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+
                     DEM, data=data, method="optimal") 
matchrun=cbind(matchrun,summary(model)$sum.matched[,c(1:2,4)])

model <- matchit(DRUGP~ALLGEMP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+
                     Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+
                     DEM, data=data, method="full") 
matchrun=cbind(matchrun,summary(model)$sum.matched[,c(1:3)])

# Table VIII balance statistics for total battle deaths
round(matchrun[-1,],2)

## Duration balance statistics (bottom panel, Table VIII)
# Select variables for the analysis and pare the dataset down
keep.vars <- c("BTTLDtotLN","BTTLDtot","resources","ongoingCONF","Duration","DRUGP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","DEM","ccode","terr","gdplLN")  
data <- na.omit(dataset[keep.vars])

# Create censored observations
data$ongoingCONF=1-data$ongoingCONF

# Run analysis
matchrun=c()
model <- matchit(DRUGP~mountain+resources+gdplLN+Weakrebel
                   +Equalstrenght+internatz+coldwar+popLN+
                     ethfrac+DEM+terr+confbord, 
                   data=data, method="nearest")
matchrun=cbind(summary(model)$sum.all[,c(1:2,4)],summary(model)$sum.matched[,c(1:2,4)])

model <- matchit(DRUGP~mountain+resources+gdplLN+Weakrebel
                   +Equalstrenght+internatz+coldwar+popLN+
                     ethfrac+DEM+terr+confbord, 
                   data=data, method="optimal")
matchrun=cbind(matchrun,summary(model)$sum.matched[,c(1:2,4)])

model <- matchit(DRUGP~mountain+resources+gdplLN+Weakrebel
                   +Equalstrenght+internatz+coldwar+popLN+
                     ethfrac+DEM+terr+confbord, 
                   data=data,  method="full")
matchrun=cbind(matchrun,summary(model)$sum.matched[,c(1:3)])

# Table VIII balance statistics for duration
round(matchrun[-1,],2)